Solitons and vortices in honeycomb defocusing photonic lattices 
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Solitons and necklaces in the first band-gap of a two-dimensional optically induced honeycomb 
defocusing photonic lattice are theoretically considered. It is shown that dipoles, soliton necklaces, 
and vortex necklaces exist and may possess regions of stable propagation through a photorefractive 
crystal. Most of the configurations disappear in bifurcations close to the upper edge of the first 
band. Solutions associated with such bifurcations are also numerically examined, and it is found 
that they are often asymmetric and more exotic. The dynamics of the relevant unstable structures 
are also examined through direct numerical simulations revealing either breathing oscillations or, in 
some cases, destruction of the original waveform. 

PACS numbers: 



I. INTRODUCTION 

In the past few years, there has been a considerable 
growth of interest in the examination of the self-trapping 
of light in photonic lattices optically induced in nonlin- 
ear photorefractive crystals, such as strontium barium 
niobate (SBN). This can be attributed to a considerable 
extent to the fact that the theoretical inception [l[ of 
the relevant phenomena was rapidly followed by the ex- 
perimental realization [H, 0, 0] ^ revealing a considerable 
wealth of new possibilities. This setting naturally per- 
mits the consideration of the competition between ef- 
fects of nonlinearity and those of diffraction, therefore 
enabling the examination of effects of periodic "poten- 
tials" on solitary waves. In this context the role of the 
effective potential is played by the ordinary polarization 
of light forming a waveguide array in which the nonlinear, 
extra-ordinarily polarized probe beam evolves. 

Numerous nonlinear waves and coherent structures 
have been elucidated and experimentally realized in this 
context. In particular, discrete dipole [H, necklace 
solitons and even stripe patterns |l3], rotary solitons [8|, 
discrete vortices |9|] or the realization of photonic qua- 
sicrystals [lo| and Anderson localization are among 
the recently reported experimental results in the field. 
These efforts illustrate the potential that this setting 
holds for the examination of localized structures that may 
be usable as carriers and conduits for data transmission 
and processing in all-optical communication schemes. In 
parallel to this more practical aspect, this framework 
remains an experimentally tunable playground where 
numerous fundamental issues of solitons and nonlinear 
waves can be explored. 

The above mentioned interplay of nonlinearity with pe- 
riodicity is important not only in the physics of optically 
induced lattices in photorefractive crystals, but also in a 
variety of other contexts in optical and atomic physics. 
These involve e.g. on the optical end, the numerous de- 
velopments on the experimental and theoretical investi- 
gation of optical waveguide arrays; see e.g. [l2|, p^3] for 
relevant reviews. In the case of atomic physics, and par- 
ticularly of Bose-Einstein condensates, the confinement 



of dilute alkali vapors in optical lattice potentials 
has offered a similarly far-reaching opportunity to exam- 
ine many fundamental phenomena involving (effective) 
nonlinearity and spatial periodicity. These include, but 
are not limited to modulational instabilities, Bloch oscil- 
lations, Landau-Zener tunneling and gap solitons among 
others; see [15] for a recent review. 

Our present study, motivated by optically induced lat- 
tices in photorefractive SBN crystals, focuses on two- 
dimensional periodic, nonlinear media with a non-square 
lattice. While most of the above studies have been ded- 
icated to square lattices, only a few have tackled the 
coherent structures possible in non- square settings; see 
e.g., as relevant examples [l6|, [13, E, Ell and ref- 

erences therein. Furthermore, the vast majority of the 
above-mentioned studies has centered around focusing 
nonlinear it ies. At least partly, this is due to technical 
limitations, as it is easier to work with voltages that are 
in the regime of focusing rather than in that of the defo- 
cusing nonlinearity (in the latter case, sufficiently large 
voltage, which is tantamount to large nonlinearity, may 
actually change the sign of the nonlinearity by inverting 
the orientation of the permanent polarization of the crys- 
tal). As a result, coherent structures in the defocusing 
regime, have only rather sparsely been examined. Such 
an experimental example is the fundamental and higher 
order gap solitons excited in the vicinity of the edge of 
the first Brillouin zone [H, Hlj . More complex gap struc- 
tures (multipoles and vortices) are only now starting to 
be explored in square lattices [22]. In parallel to these 
experimental developments, a theoretical framework is 
starting to emerge to address such multipole and vor- 
tex structures in square lattices with cubic nonlinearities 
p3l , [2^ , whose qualitative predictions can however be ex- 
tended to non-square settings and the main ones among 
which will also be compared to the results presented be- 
low. Our main focus in the present work is on employing 
a continuum model to examine the waveforms present in 
a context involving a triangular lattice {honeycomb ) po- 
tential and a saturable defocusing nonlinearity associated 
with appropriate optically induced lattices in SBN crys- 
tals. In particular, we study in detail multipole (dipole 
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and hexapole) solitons in such lattices induced with a 
self-defocusing nonhnearity. 

We numericahy analyze both the existence and the sta- 
bility of these structures and follow their dynamics, in the 
cases where we find them to be unstable. We also quali- 
tatively compare our findings with the roadmap provided 
by the discrete model [23]. 

Our presentation is structured as follows. In section 
II, we present our theoretical model setup. Dipole solu- 
tions with the two excited sites in adjacent wells of the 
periodic potential (nearest-neighbor dipoles) are studied 
in section III. Subsequently, we do the same for next- 
nearest-neighbor dipoles (excited in two diagonal sites, 
separated by one lattice site) in section IV and opposite 
dipoles on either end of the hexagonal configuration (ie. 
the excited sites are separated by two empty wells) in 
section V. Section VI addresses the case of more com- 
plex structures such as hexapoles (all six sites from one 
period of the potential) and vortices. Finally, in section 
VII, we summarize our findings, posing some interesting 
questions for future study. 



II. SETUP 

We use the standard partial differential equation for 
the amplitude of the electric field U [19,, llO, [H, H , in 
the following form: 



iU, = [L + N{x,\U\^)]U, 



7V(x,|C/|2) = 



En 



./(X) + |C/|2' 



(1) 



(2) 



where L = DV^ and is the two-dimensional Lapla- 
cian, U is the slowly varying amplitude of the probe 
beam, and 



/(x) =/o|e' 



Jfebix 



(3) 



is the optical lattice intensity function formed by three 
laser beams with bi = (1,0), b2 = {.— \^—^)^ and 

ba = (— ^,-^). Here /q is the lattice peak inten- 
sity, z is the propagation distance and x = {x^y) are 
transverse distances (normalized to = 1 mm and 
Xs = Us = 1/im), £^0 is proportional to the applied DC 
field voltage, D = Zs\/ {^"KneXsHs) is the diffraction co- 
efficient, A is the wavelength of the laser in a vacuum, d 
is the period in the x direction with k = 47r/(3(i) (period 
in the y direction is a/S^), and Ue is the refractive index 
along the extraordinary axis. We choose the lattice in- 
tensity /o = 0.6. A plot of the optical lattice is shown 
in Fig. [1] for illustrative purposes regarding the location 
where our localized pulses will be "inserted" . In addition, 
we choose other physical parameters consistently with a 
typical experimentally accessible setting [22[ as 



d = 30/im, A = 532 nm, rie = 2.35, Eq = 8. 

The non-dimensional value D = 18.01, and we note 
that this dispersion coefficient is equivalent to rescaling 
space by a factor ^/D as e.g. in [27]. 

The numerical simulations are performed in a rectan- 
gular domain corresponding to the periodicity of the lat- 
tice using a rectangular spatial mesh with Ax ^0.75 and 
Ay ^ 0.86 and domain size 4d X 3 V3(i, i.e. 160 x 180 grid 
points. See Fig. [1] for a schematic of the spatial configu- 
rations. 

Regarding the typical dynamics of a soliton when it is 
unstable, we simulate the z-dependent evolution using a 
Runge-Kutta fourth-order scheme with a step Az = 0.01. 
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FIG. 1: (Color online) A spatial (x-y) contour plot of the or- 
dinary polarization standing wave [lattice beam in Eq. (|3])]. 
In this context, the light intensity maxima correspond to the 
minima of the resulting refractive index lattice (i.e., honey- 
comb lattice), as opposed to the focusing nonlinearity lattice 
field, where they correspond to the maxima (i.e., triangular 
lattice). Points A, C, D, E, and F are used for naming 
the various configurations. A is a "nearest-neighbor" mini- 
mum of B and F, a "next-nearest-neighbor" of C and E, and 
an "opposite" of D (with respect to the local maximum of 
the lattice). Because of the symmetry of the setup, this is 
a complete characterization of dipole configurations. We will 
refer to the configurations with the names given above. 

Assuming a stationary state u{x^ y) exists, and letting 
the propagation constant ji represent the (nonlinear) real 
eigenvalue of the operator of the right-hand-side of Eq. 
([1]), then the corresponding eigenvector u{x,y) is a fixed 
point of 



[/i-L-N{^,\u\^)]u = 0, 



(4) 



The localized states it of (|4]) were obtained using the 
Newton- GMRES fixed point solver nsoli from [28] and 
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a pseudo arc-length continuation [29] was used to follow 
each branch and locate the bifurcations which occur at 
the edge of the first band. Since the parameter of interest 
is /i, diagnostics are plotted against ja. 

We restrict /i to those values within the first spectral 
gap of the linear eigenvalue problem, 



[/i-L- Ar(x,0)]^i = 0. 



(5) 



Values of the propagation constant /i within this for- 
bidden gap in the spectrum of the linearized problem will 
correspond to exponentially localized in space, so-called 
gap-soliton, states of the original nonlinear partial dif- 
ferential equation. Using a standard eigenvalue solver 
package implemented through MATLAB, we identify the 
spectral gap for our given parameters and gridsize to be 
3.62 < /i < 4.94. 

The square root of the optical power (or, mathemat- 
ically, the norm) of the localized waves is defined as 
follows: 



P = 



\U\^dxdy 



1/2 



(6) 



Introducing a linearization around an exact stationary 
solution li, and expanding the leading order perturba- 
tion into a eigenfunctions and eigenvalues, we obtain the 
following Bogoliubov system 



zA + /i — L 
iX — fi -\- L 



d{NU) 
dU ' 
d{NUy 



d{NU) 

U 7^TT—\uU = 0, 



d{NUy 
" dU 



\uu = 0.(7) 



We solve the above linear eigenvalue problem using MAT- 
LAB's standard eigenvalue solver package. The symplec- 
tic nature of the resulting eigenvalue equations guaran- 
tees that the relevant eigenvalues should come in quar- 
tets, hence an instability is present whenever the solution 
of the above linearization problem of Eqs. ([7j) possesses 
an eigenvalue with a non-zero real part. 

We now briefly discuss the principal stability conclu- 
sions, for the defocusing case of which we should ex- 
pect to still be valid in the present configuration. Nearest 
neighbor excitations in the defocusing case correspond 
to nearest neighbor excitations in the focusing case, but 
with an additional tt phase in the relative phase of the 
sites added by the so-called staggering transformation 
[23]. Therefore, the in-phase nearest neighbor config- 
uration in the defocusing case corresponds to an out- 
of-phase such configuration in the focusing case (and 
should thus be stable) [2?]. On the other hand, next 
nearest neighbor out-of-phase defocusing configurations 
would correspond to next nearest neighbor out-of-phase 
focusing configurations and should also be stable (at least 
in some parameter regimes). By the same token, out-of- 
phase nearest neighbor, and in-phase next nearest neigh- 
bor structures should be unstable. These considerations 



also indicate that in-phase opposite dipoles should be 
stable, while out-of-phase such dipoles should always be 
unstable. Finally, vortex-like structures and in-phase 
hexapoles should be stable as well. Notice, however, that 
as discussed in [23] the multipole structures characterized 
as potentially stable above will, in fact, typically possess 
imaginary eigenvalues of negative Krein signature (see 
e.g. [30] and references therein). These may lead to os- 
cillatory instabilities through complex quartets of eigen- 
values. These arise by means of Hamiltonian-Hopf bifur- 
cations [3l| emerging from collisions with eigenvalues of 
opposite (i.e., positive) Krein signature. These conclu- 
sions will be discussed in connections with our detailed 
numerical results in what follows. 



III. NEAREST NEIGHBOR DIPOLE SOLITONS 

In this section, we report dipole solitons where the two 
lobes of the wave are located in two nearest neighbor 
(N) lattice sites in the 2D triangular potential shown in 
Fig. [H The lobes can have the same phase or tt phase 
difference so we define them as in-phase (IP) dipoles and 
out-of-phase (OP) dipoles, respectively. 



A. In-Phase Nearest Neighbor Dipole Solitons 

We have found IP dipoles in adjacent wells for values of 
the propagation constant fi throughout the entire Bragg 
reflection gap for a given ^o- We found that the solitons 
exist for /i between 3.62 and 4.94, and that the intensity 
of the dipoles cannot be arbitrary low, a result similar to 
the observed results of the focusing and defocusing cases 
for square lattices [H, IH, [27|] . The relevant findings are 
summarized in Fig. [2l 

The top left panel of Fig. [2] shows the stability of the 
dipoles against the propagation constant /i, by illustrat- 
ing the maximal growth rate (maximum real part of all 
eigenvalues A) of perturbations. When max(Re(A)) = 0, 
this implies stability of the configuration, while the con- 
figuration is unstable if max(Re(A)) 7^ in this Hamilto- 
nian system. We found that this type of dipoles may be 
stable for windows throughout the first Bragg gap, as pre- 
dicted above, although it is possible for small oscillatory 
Hopf instabilities to arise due to opposite signature eigen- 
value collisions. The dipole configuration disappears in a 
saddle-node bifurcation at the edge of the first spectral 
band, depicted in the top panels of Fig. [21 as /i ^ 4.94, 
and a real pair of eigenvalues emerges. At this point, 
the configuration collides with a configuration shown at 
the bottom panel of Fig. [2] in which the adjacent well 
next to one of the populated ones becomes excited out- 
of-phase with the others. Consistent with our theoretical 
expectation from its having an out-of-phase set of near- 
est neighbors, the latter configuration always has a real 
pair of eigenvalues A. 
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FIG. 3: The typical time- dependent dynamics of an unstable 
configuration along the upper (dashed line) branch of the ex- 
istence curve presented in the top panels of Fig [2] Depicted 
here is the isosurface of height 0.15 of the dynamics of the of 
the intensity, of the configuration shown in the bottom 

panel of Fig. O 



FIG. 2: (Color online) The top left panel shows the stabil- 
ity of the dipoles against the propagation constant /i. It 
is stable when the spectra is purely imaginary (i.e., when 
max(Re(A)) = 0). The top right panel depicts the power 
of the dipoles against the propagation constant. In each of 
these images the solution branch is denoted by a solid line. 
The branch with which the dipole collides and terminates in 
a saddle-node bifurcation is shown by a dashed line. The 
shaded areas in both of these panels represent the bands of 
linear spectrum (|5}. The middle left and right panels show 
the profile u of the dipole at /x = 4 and the corresponding 
complex spectral plane (Re(A), Im(A)) of A = Re(A) +iIm(A). 
Finally, the bottom panels show the same features for the 
unstable saddle solution corresponding to the dashed line. 



The middle left and right panels show the profile u 
of the in-phase nearest (IPN) neighbor dipole at /i = 4 
and the corresponding spectrum of linearization eigen- 
values A = Xr -\- iXi in the complex plane (Ar,Af), re- 
spectively. The corresponding profile and spectral plane 
for the saddle branch (that eventually collides with the 
IPN solution) at /i = 4 is shown in the bottom left and 
right panel, respectively, of the same figure, illustrating 
the exponential instability of the latter. 

We have simulated the dynamics of the solitary waves 
when they are unstable. The dipoles are perturbed by 
a random noise with maximum intensity 2 x 10~^. It 
is interesting to note that an unstable IPN dipole turns 
out to be quite robust, even though it experiences only 
an oscillatory instability. It is remarkable that up to 
z = 200 we did not see any signifant change in the con- 
figuration. Therefore, we do not depict our evolution 
simulation here; we simply note that this is consonant 
with the very weak growth rate of the relevant oscilla- 
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FIG. 4: (Color online) The top panels correspond to the same 
panels of Fig. [2] but for OPN dipole solitons. The bottom 
panels show the profile u and the corresponding spectral plane 
of the dipoles at /x = 4. 



tory instability. 

For the solution branch shown in the bottom panel of 
Fig. [21 we present its dynamics in Fig. [3l We found that 
the instability is strong as predicted above such that even 
after a relatively short propagation distance, the insta- 
bility already sets in and leads to recurrent oscillations 
(for the remainder of our dynamical evolution horizon) 
between a dipole, two-site state and a three-excited-site 
state; see Fig. [3l 
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FIG. 5: The typical dynamical evolution of an unstable out- 
of-phase nearest neighbor configuration from the family pre- 
sented in FiglH Depicted is the isosurface at half the maxi- 
mum of the intial intensity amplitude. Notice that the OPN 
appears to oscillate between two sites and one for the prop- 
agation constant /x = 4 (center) of the solution presented 
in the bottom panels of Fig. IH as does it for smaller values 
(/i = 3.6, left), although for a larger value of /x = 4.6, right, 
the solution essentially transforms (due to the instability) into 
a single site mode. We hypothesize that the stronger effect 
of the nonlinearity on those solutions with smaller values of 
/i decreases the size of linear regime, causing these more lin- 
early unstable solutions to actually appear more stable in the 
full nonlinear dynamical evolution, although it is away from 
the linear regime that the instability is manifested merely as 
oscillations. 



B. Out of Phase Nearest Neighbor Dipole Solitons 

We have also found OP dipoles arranged in nearest- 
neighboring lattice wells which we refer to as OPN. We 
summarize our findings in Fig. [H where one can see that 
the solitons exist in the whole entire region of propaga- 
tion constant /i in the first Bragg gap, /i G (3.62,4.94). 
This smooth transition indicates that the OPN dipole 
solitons emerge out of the Bloch band waves; see e.g. (32j 
and [s^ for a relevant discussion of the ID and of the 2D 
problem respectively, in the case of the cubic nonlinearity. 
The OPN dipoles are unstable due to a real eigenvalue 
pair, as expected from our above theoretical predictions. 

As the branch merges with the band edge, we observe 
an interesting feature, namely that the configuration be- 
gins to resemble a hexapole with a tt phase difference 
between each well. This can be an indication that these 
structures bifurcate out of the Bloch band from the same 
bifurcation point. We elaborate this further in our dis- 
cussion at the end of section VL 

In Fig. [5] we present the unstable dynamics of OPN 
dipole solitons perturbed by similar random noise per- 
turbation as in Fig. [3l We display here three solutions 
for a range of chemical potentials to illustrate that the 
dynamical evolution of linearly unstable states is appar- 
ently correlated to the power of the solution. This type 



of dipoles is typically more unstable than its IP counter- 
part, as is illustrated in the figure. In particular, in all 
three examples of unstable evolution given the instability 
already starts to manifest itself, around z — 20 However, 
for small values of fi (large power) the OPN continues 
oscillating between a single site structure and a two site 
structure for the (longer) evolution distances investigated 
in this illustrative case, while for large enough /i (small 
enough power), one of the sites decays and the power is 
concentrated on a single site. 



IV. NEXT NEAREST NEIGHBOR DIPOLE 
SOLITONS 

We have also obtained dipole solutions that are not 
oriented along the two nearest-neighboring lattice wells, 
but rather where the two humps of the structure are lo- 
cated at two next-nearest-neighboring lattice sites. These 
humps can once again have the same phase or a tt phase 
difference between them. We will again use the corre- 
sponding IP and OP designations for these next nearest 
neighbor waveforms. 



A. In Phase, Next Nearest Neighbor Dipole 
Solitons 

The in-phase next-nearest (IPNN) neighbor solitons 
exist only up to a marginal distance from the second 
band. The stability and power of these dipoles are shown 
in Fig. [6l The stability is again consistent with the the- 
oretical discussion of Section II. In particular, the IPNN 
configuration always possesses a real eigenvalue pair; fur- 
thermore, the corresponding unstable "saddle" structure 
with which it collides and terminates through a saddle- 
node bifurcation has an additional such eigenvalue pair 
(two real eigenvalue pairs in total for the solution branch 
indicated by dashed line in Fig. [6|). 

We have simulated also the dynamics of the unstable 
IPNN. Yet, we do not present our simulation here as the 
typical evolution of this configuration is quite in resem- 
blance to the dynamics of an unstable OPN (see Fig. [5]) in 
the fact that the configuration recurrently oscillates be- 
tween a two-soliton state and a one-soliton state. Such 
an oscillation persists even up to z = 200. 

In Fig. [71 we present the dynamical evolution of the bi- 
furcating solution shown in the bottom panel of Fig.[6]un- 
der similar random noise perturbation as above. One can 
note similarities in the typical evolution of this configu- 
ration and the evolution of the bifurcating IPN solution 
shown in Fig. [3l one of which is the recurrent oscillation 
between a pattern with three pulses and one with just 
two peaks. 
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FIG. 6: (Color online) The top panels correspond to the same 
diagnostics as in Fig.[2]but for IPNN dipole solitons. The mid- 
dle panels show the profile U and the corresponding spectral 
plane of the IPNN dipole at /i = 4.1, while the bottom row 
shows the same images for the solution branch corresponding 
to the dashed line in the top panel, shown at the same value 
of /i. 




FIG. 7: The same figure as Fig. O but for the solution pre- 
sented in the bottom panel of Fig.[6l Depicted is the isosurface 
of height 0.05. 



B. Out of Phase Next Nearest Neighbor Dipole 
Solitons 

We have also obtained out-of-phase, next-nearest 
(OPNN) neighbor dipole solitons. A typical profile of 



FIG. 8: (Color online) The top panels depict the largest real 
part of the critical eigenvalue, as well as the power of the 
OPNN dipole solitons. The middle panels show the profile 
u and the corresponding spectra in the complex plane of the 
dipole at /i = 3.9, and the bottom is the unstable saddle 
configuration at the same value of /i, where one of the sites 
has merged with a neighbor out of phase and become an OPN, 
accounting for the real eigenvalues. 




FIG. 9: The same figure as Fig. O but for the solution pre- 
sented in the bottom panel of Fig. [8] Depicted is the isosurface 
of height 0.05. 

this family of solutions for /i = 3.9 is shown in Fig. [HI 
The power diagram of these solitons is presented in the 
top panel of Fig. [51 Typically these structures are sta- 
ble (as indicated again by the comparison with the the- 
oretical discussion and by the numerical results shown 
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in the middle right panel of Fig. [8] ), suffering only win- 
dows of oscillatory instability due to the presence of a 
single eigenvalue with negative signature and its colli- 
sion with the spectral bands In fact, we have found that 
a consistent stability range for Eq = S exists between 
4.5 </i< 4.85. 

For this solution we also observe that, similarly to the 
IPNN dipoles, the solution disappears at non-zero inten- 
sity because of the collision of this dipole with another 
(three-site) configuration shown in the bottom panels of 
Fig. [8] in a saddle- node bifurcation. It is relevant to note 
that the point of the bifurcation is very close to the edge 
of the Bloch band, i.e., to /j. ^ 4.94. 

The dynamics of the OPNN dipole do not manifest 
their very weak oscillatory instability for the evolution 
distances considered herein. On the other hand, the dy- 
namics of the instability of the three-site solution (with 
which the OPNN branch collides in the saddle-node bi- 
furcation) can be seen in Fig. [9l More specifically, the 
instability manifests itself in the form of interactions be- 
tween the closest out-of-phase pair of solitons (leading 
to recurrent oscillations between a three-peak and a two- 
peak state). Notice that the third peak is almost not 
affected by these interactions. 

V. OPPOSITE DIPOLE SOLITONS 

We now address opposite (O) dipole solitons residing 
at the two sites along a diameter of a local maximum of 
the lattice. This is the final type of dipole configuration 
for a symmetric triangular lattice, exhausting the pos- 
sibilities up to phase and rotational invariances. Again, 
we partition our considerations into in-phase and out-of- 
phase cases. 

A. In Phase Opposite Dipole Solitons 

We have found in-phase opposite (IPO) solitons 
throughout the first gap in the linear spectrum. Our 
numerical findings are presented in Fig. [TOl 

Again, the solution branch is largely stable with small 
windows of Hopf quartets and again a saddle node bifur- 
cation occurs as the branch approaches the first spectral 
band. Also, interestingly, the configuration with which 
this branch collides when it disappears resembles an OPN 
(or two pairs of OPNs- see the third and fourth row of 
the figure). The latter branches are naturally unstable 
due to one (or more) real pair of eigenvalues. 

The dynamics of one of the bifurcating solutions, i.e. 
the configuration with a single OPN structure, is pre- 
sented in Fig. [TTl where one can see that, as usual, only 
the pair of out-of-phase nearest neighbor dipole interacts, 
while the other soliton is almost uninfiuenced. 

Using the same reasoning, one can deduce as well that 
the dynamics of the other bifurcating solution, presented 
in the bottom panel of Fig. [TOl will be similar, except 
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FIG. 10: (Color online) The top panels depict the largest 
real part of the critical eigenvalue, as well as the power and 
the peak intensity of the IPO dipole solitons. The panels 
in the second row show the profile u and the corresponding 
spectra in the complex plane of the dipole at /x = 4.5, the 
third row shows the same images at the same value of fi for 
the middle branch (dashed line) of the bifurcation diagram 
and the bottom row is a solution along the top branch (dash- 
dotted line) at the same value. 



the fact that now there are two pairs of OPN interacting 
among themselves. 



B. Out of Phase Opposite Solitons 

Lastly, as regards dipoles, we consider the out of phase 
opposite (OPO) dipole. The first interesting characteris- 
tic of the OPO is its strong instability stemming from a 
real pair of eigenvalues, seen in the top left and middle 
rows of Fig. [121 Once again the direct instability of 
this mode follows from our theoretical considerations of 
Section II. On the other hand, the figure also reveals an 
interesting bifurcation structure in this case. The branch 
actually merges with a hexapole made of three copies of 
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FIG. 11: The same figure as Fig. O but for the solution pre- 
sented in the middle panel of Fig. IIQI Depicted is the isosur- 
face of height 0.05. 
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FIG. 12: (Color online) The top panels depict the largest real 
part and the power of the OPO dipole solitons. The middle 
panels again show the profile u and spectra at /x = 4, and the 
bottom is the more unstable saddle configuration, consisting 
this time of a hexapole configuration constructed out of three 
such OPOs. 



itself close to the band, when solutions start becoming 
extended. This hexapole then intersects with the linear 
spectrum shortly thereafter and the solution transforms 
itself into a fully extended "checkerboard" -like configu- 
ration of all adjacent wells excited out-of-phase. As can 
be seen in the top left and the bottom right of Fig. [121 
the hexapole configuration is significantly more unstable, 



FIG. 13: The top panel is the same as Fig.[3l but for the solu- 
tion presented in the middle panel of Fig. 1121 with isosurface 
of height 0.1. On the other hand, the bottom panel illustrates 
again (c.f. Fig. [5]) that the linear stability analysis is more 
predictive of the nonlinear dynamics for the solution with the 
larger value of /i = 4.88 (and accordingly smaller amplitude) 
close to the intersection with the extended OP quadrupole 
branch (the isosurface is taken at half the maximum of the 
initial intensity amplitude). The growth rates for each so- 
lution are comparable, while the dynamical evolutions differ 
drastically. 



possessing five real eigenvalue pairs. 

We have numerically monitored the full evolution to 
observe the dynamics of the unstable OPO dipoles. It is 
interesting to note that even though the state has a pair 
of real eigenvalues, our simulation reveals that the insta- 
bility is barely detectable for the state depicted in the 
middle rows of Fig. [121 presumably because of the spa- 
tial separation of the populated sites (top row of Figure 
[T3|) : the solution oscillations are very mild (and almost 
indetectable) between similar structures with mass con- 
centrated in one site or another. On the other hand, for 
significantly smaller power (larger /i) as seen in the bot- 
tom panel of Figure [131 oiie site decays fairly rapidly 
and a robust single site remains. 

Regarding the bifurcating solution, which is an out-of- 
phase hexapole, we will explore it as well as the other 
hexapole configurations in more detail in the following 
section. 
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FIG. 14: The same dynamical evolution figure as Fig. El but 
for the out-of-phase hexapole depicted in the bottom panel of 
Fig. [12] Depicted is the isosurface of height 0.1. 
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FIG. 15: (Color online) The top panels depict the largest real 
part and the power of the IP hexapole solitons. The middle 
panels again show the profile u and spectra at /i = 4, and 
the bottom is the more unstable saddle configuration, which 
features our expected OPN sidekick. 



VI. HEXAPOLE SOLITONS AND VORTEX 
NECKLACES 

First, we consider the out-of-phase hexapole. The ex- 
istence and the stability of this configuration has been 
described in the preceding section. As the state has mul- 
tiple pairs of real eigenvalues, it is natural to expect that 



FIG. 16: The same figure as Fig. O but for the in-phase 
hexapole depicted in the bottom panel of Fig. 1151 Depicted 
is the isosurface of height 0.3. 



it should be prone to break up under the instability's dy- 
namical evolution. A typical example of such a numerical 
experiment is presented in Fig. [Ml 

We found IP hexapole configurations as well, which, 
in accordance with our considerations in Section II, turn 
out to chiefly be stable within the first gap, although 
they may possess weak oscillatory instability inducing 
eigenvalue quartets. 

This configuration also suffers a saddle-node bifurca- 
tion with an OPN-type pair emanating off of one of its 
lobes, when a neighboring well becomes populated out 
of phase near the first band. The latter configuration is 
unstable always possessing a real eigenvalue pair in its 
linearization spectrum. We note in passing that this is 
among any of the six equivalent symmetric versions of 
this configuration. 

As for the dynamics of the instability, the solution 
along the main lower branch is quite robust to strong per- 
turbation. Even though the solution suffers from an os- 
cillatory instability, a random perturbation with a max- 
imum intensity almost 10~^ cannot lead to a breakup of 
the configuration until propagation distances of the order 
of z = 200. On the other hand, the oscillatory dynam- 
ics leading to the break up of the configuration of the 
bottom panel of Fig. [15] is shown in Fig. [161 

Finally, we investigate the complex-valued hexapole 
configuration for which each lobe has the same modu- 
lus and their phase increases counterclockwise in phase 
increments of 7r/6, yielding a vortex- necklace configura- 
tion. This configuration turns out to be stable for the 
most part within the first gap as well, with minor Hamil- 
tonian Hopf-bifurcation induced oscillatory instabilities. 
We also found that this solution undergoes a saddle-node 
bifurcation near the first band, in which it collides with 
a waveform with two pairs of OPNs. The stability of the 
latter configuration in the presence of these additional 
OPN dipoles is consistent with that of their real coun- 
terparts from the previous sections, each appearing to 
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FIG. 17: (Color online) The top panels depict the largest real part and the power of the vortex necklace hexapole configuration. 
The second row depicts the modulus of the solution and corresponding spectrum when /a = 4.6. The third row illustrates 
the real and imaginary components of the field and the phase (from left to right). The fourth and fifth rows show the same 
properties as the second and third but for the unstable eight-site configuration of the dashed line in the top panels. 
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FIG. 18: The same figure as Fig. O but for tiie vortex neck- 
lace witii eigiit lobes depicted in the bottom panel of Fig. [T71 
Depicted is the isosurface of height 0.1. 

contribute one real pair, rendering the relevant configu- 
ration quite unstable. 

Similar to the case of in-phase hexapoles, even though 
vortex necklaces may be unstable, they are quite robust 
to perturbation, given the weak nature of the relevant 
oscillatory instabilities. We therefore only depict the dy- 
namics of the solutions which have eight lobes as shown 
in the fourth and fifth row panels of Fig. [T71 The typical 
evolution of this state is shown in Fig. [181 showcasing 
the oscillatory breakup of this structure into one with a 
smaller number of lobes. 



VII. CONCLUSIONS 

In this communication, we examined in detail theo- 
retically and numerically the existence, stability and dy- 
namics of multipole lattice solitons excited with a sat- 
urable defocusing photorefractive nonlinearity in a tri- 
angular geometry. We have obtained a wide array of 
relevant structures, including different types of dipoles 
and hexapoles, as well as vortices. For the dipole con- 
figurations we examined the different possible phase con- 
figurations (in phase, and out phase profiles), as well as 
cases where the excited sites are separated by 0, 1, or 
2 intermediate lattices sites. For hexapoles, we exam- 
ined in phase and out of phase, and we also studied the 
monotonic increasing phase of a discrete vortex necklace. 

We have found good agreement with the general guide- 



lines, explained in section II, stemming from the theoret- 
ical analysis of the discrete model. This intuition led to 
the illustration of a wide variety of potentially stable so- 
lutions (although they may incur oscillatory instabilities) 
such as the in phase, nearest neighbor dipole, the out of 
phase, next nearest neighbor dipole, and the in phase 
opposite dipole. We have also identified those solutions 
including e.g., the out of phase nearest neighbor, in phase 
next-nearest neighbor, and out of phase opposite dipoles 
which are typically unstable due to exponential instabili- 
ties and real eigenvalues. By the same considerations, the 
in-phase hexapole was proposed and was indeed found to 
be typically stable, while the out-of-phase one was pre- 
dicted and observed to be quite unstable, due to multiple 
real eigenvalue pairs. Finally, we have seen that the dis- 
crete vortex structure is also potentially stable. 

Furthermore, we have also identified an interesting set 
of bifurcations that are associated with the paramet- 
ric continuation and termination of some of the above 
branches. The dynamical instabilities encountered in the 
present work have been monitored through direct inte- 
gration of the relevant dynamical equation. The result of 
evolution in every case involved oscillations between the 
original configuration and one with fewer sites which is 
more stable, such as a single site solitary wave, and some- 
times degeneration to such a configuration. Solutions 
with smaller power tend to decay into a single site soli- 
tary wave for certain solutions investigated, while those 
with larger power tend to oscillate. This connection is 
beyond the scope of the present work, but is currently 
being investigated further. 

Since the framework of defocusing equations has been 
studied far less extensively than their focusing counter- 
parts, it would be particularly interesting to extend the 
present considerations to other structures. Perhaps the 
most interesting example would be the study of multiple 
charge vortices in this context which would be an inter- 
esting endeavor both from a theoretical, as well as from 
an experimental point of view. Such studies are currently 
in progress and will be reported in future publications. 
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